function [rmBlocks,rateMatching] = WiSiL_Encoder_rateMatching(...
                                        encodedBlocks,transportBlockSize,wordSize,layers)

% This function should be executed after the turbo coding and before the
% block concatenation. 
%
% Inputs:
% encodedBlocks: A cell array 1 x C containing blocks of encoded bits. Each
% element of the array is a 3 x D matrix and represents the output of the
% turbo encoder. The first row contains the systematic bits, the second row
% contains the 1st encoder parity bits and the third row contains the
% 2nd encoder parity bits.
%
% transportBlockSize: it is the size of the output vector of the function
% "WiSiL_Tx_channelCoding. It is used to determine the size of the output
% vector E.
%
% wordSize: it is the size in bits of the modulated symbols
%
% layers: number of layers related with MIMO transmissions
%
% Outputs:
% rmBlocks: A cell array 1 x C. Each cell element is a 4 X E matrix. Each
% line of the matrix is an output of the rate matching process depending on
% the redundancy version rvdix.
%
% rateMatching: a structure containing important information about the rate
% matching process to be used at the decoder.
%
% The function performs the following tasks:interleaving,bit collection and
% bit selection according to the specification "3GPP TS 36.212 V12.1.0".
% The reader is referred to this document for a more detailed description 
% of what the code does.
%
% History:
%
%   Lucas Sousa e Silva 08/2014 - created


%-----------------------Sub Block interleaver-----------------------------
rmBlocks = cell(7,1);
P = [1,17,9,25,5,21,13,29,3,19,11,27,7,23,15,31,2,18,10,26,6,22,14,30,4,...
    20,12,28,8,24,16,32];
    rvdix = [0 1 2 3];                                                      % Redundancy version for HARQ
C_subBlock = 32;                                                            % Number of columns.
D = zeros(1,length(encodedBlocks));                                         % Number of row elements on each encoded Block
N_cb = zeros(1,length(encodedBlocks));                                      % Size of the soft buffer
k0 = zeros(length(encodedBlocks),length(rvdix));                            % offset of the output buffer (HARQ)

for r = 1:length(encodedBlocks)
    
    D(r) = size(encodedBlocks{r,1},2);                                      % Length of the bitstream.
    R_subBlock = ceil(D(r)/C_subBlock);                                     % Number of rows.
    N_D = R_subBlock*C_subBlock - D(r);                                     % Number of dummy bits.
    v = zeros(3,C_subBlock*R_subBlock);                                     % Initialization of matrix v.
    
    for t = 1:3
        y = [nan(1,N_D) encodedBlocks{r,1}(t,:)];                           % Dummy bits(nan)are inserted at the beggining of the bitstream.
        if (t < 3)
            Y = vec2mat(y,C_subBlock);                                      % Writes the bitstream into a matrix R_subBlock X C_subBlock.
            Y_p = Y(:,P);                                                   % Permutes the columns according to vector P.
            v(t,:) = Y_p(:)';                                               % Output of the interleaver.
        else
            for k = 0:length(y)-1
                v(t,k+1) = y(mod((P(floor(k/R_subBlock)+1)+ C_subBlock*mod(k,R_subBlock)+1),length(y))+1);
            end
        end
    end
%----------------- Bit collection in circular Buffer-----------------------
    w = zeros(1,3*size(v,2));
    w(1:size(v,2)) = v(1,:);
    w(size(v,2)+1:2:end) = v(2,:);
    w(size(v,2)+2:2:end) = v(3,:);

% ---------------------------Bit selection---------------------------------
% K_soft,K_mimo,M_dlharq,Q_m,N_l are functions of the simulation parameters
% G is the length of the input vector of 'WiSiL_Tx_channelCoding'.
 
    K_w = length(w);
    C = length(encodedBlocks);  
    N_soft = 250368;                                                        %category 1
    if (N_soft == 35982720)
        K_c = 5;
    elseif (N_soft == 3654144 )
        K_c = 2;
    else
        K_c = 1;
    end
    K_mimo = 1;                                                             %transmission mode 1
    M_dlharq = 8;                                                           %FDD configuration
    M_limit = 8;
    N_ir = floor(N_soft/K_c*K_mimo*min([M_dlharq M_limit]));
    N_cb(r) = min([floor(N_ir/C) K_w]);   
    G = transportBlockSize;                                                     
    Q_m = wordSize;                                                         % word size according to the modulation order
    N_l = layers;                                                           % number of layers
    G_line = G/(N_l*Q_m);
    gamma = mod(G_line,C);
    
    if (r <= C - gamma - 1 +1)
        E(r) = N_l*Q_m*floor(G_line/C);
    else
        E(r) = N_l*Q_m*ceil(G_line/C);
    end

    k0(r,:) = R_subBlock*(2*ceil(N_cb(r)/8*R_subBlock)*rvdix + 2)+1;       % +1 at the end because in MATLAB index begins in 1 not 0
    e = zeros(4,E(r));
    
    for t = 1:length(rvdix)
        k = 0;
        u = 0;
        while (k < E(r))
            if (isnan(w(mod((k0(r,t)+u),N_cb(r))+1))==0)                   % +1 at the end because in MATLAB index begins in 1 not 0
                e(t,k+1) = w(mod((k0(r,t)+u),N_cb(r))+1);
                k = k + 1;
            end
            u = u + 1;
        end
    end
    rmBlocks{r} = e;
    
end

rateMatching = struct('D',D,'N_cb',N_cb,'k0',k0,'E',E);                    % Information that will be used at the decoder






